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2; Abstract 

, This paper is a slightly modified version of the introductory part of a doc- 

O toral dissertation also containing the articles hep-ph/0303076, hep-ph/0409328 and 

hep-ph/0409058. The paper focuses on the calculation of particle production in a 
relativistic heavy ion collision using the McLcrran-Vcnugopalan model. The main 
part of the paper summarizes the background of these numerical calculations. First 
we relate this calculation of the initial stage af a heavy ion collision to our under- 
standing of the whole collision process. Then we discuss the saturation physics of the 
■ small x wavefunction of a hadron or a nucleus. The classical field model of Kovner, 

McLerran and Weigert is then introduced before moving to discuss the numerical 
algorithms used to compute gluon and quark pair production in this model. Finally 
we shortly review the results on gluon and quark-antiquark production obtained in 
the three articles mentioned above. 
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Chapter 1 

Introduction 



Quantum Chromodynamics (QCD) has long since been firmly established as the correct 
theory of strong interactions; the force binding quarks into hadrons. The traditional phe- 
nomenological applications of weak coupling calculations have been in processes where 
there is large high (transverse) momentum scale given by one of the scattering particles, 
i.e. the virtual photon in deep inelastic scattering or a jet in pp collisions. This hard 
scale ensures that the relevant value of the running coupling is small enough and higher 
twist corrections are suppressed. The cross section can then be factorized into univer- 
sal parton distribution functions and perturbatively calculable partonic cross sections. 
"Bulk" observables, such as total multiplicities, have been considered so dependent on 
strong coupling physics that they can, even in theory, only be calculated from first 
principles by lattice calculations. 

The Relativistic Heavy Ion Collider (RHIC) at Brookhaven (see e.g. (U EE d for 
reviews of experimental results) has been in operation since 2000. It collides heavy ions 
and also lighter nuclei, deuterons and protons at different center of mass energies up 
to \/s = 200A GeV. One could estimate that the production of partons with trans- 
verse momentum of, say, 1 GeV at central rapidities probes the nuclear wavefunction at 
Bjorken x values of x ~ 10~ 2 . At the LHC, hopefully starting it operations in 2007 with 

= 5500 A GeV, the corresponding estimate would be x ~ 2 • 10 -4 . The mean trans- 
verse momentum of produced particles, mostly pions, at RHIC is ~ 0.5 GeV and, as we 
will discuss in Sec. l2.3( it is not unreasonable to assume the mean transverse momentum 
of the partons produced in the initial stage of the collision to be ~ 1 — 2 GeV. This 
means that the total multiplicity and transverse energy at RHIC and LHC are quantities 
that are both dominated by small x physics and involve large enough momenta to justify 
using weak coupling calculations. 

Because the constituents of a nucleus are Lorentz-contracted to the same transverse 
plane when boosted to high energy, one expects effects arising from the high density 
of "wee" partons to appear at higher values of x for larger A. One could argue that 
if nonlinear effects due to the high density of partons in the proton wavefunction start 
to be important at x < xq (for fixed Q 2 ), the corresponding effects for nuclei should 
be seen for x < Axq 1 . This estimate would mean that for the purpose of observing 
nonlinear effects due to high parton density at e.g. transverse momenta of 1 GeV RHIC 
would correspond to x ~ 3 • 10~ 5 and the LHC to x ~ 10~ 6 in deep inelastic scattering 

lr The scaling in the saturation model would be A 1 ^ 3X with A « 0.28. See e.g. |Sj for a discussion of 
the ^4-scaling in the saturation model for small x deep inelastic scattering. 
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on protons. Comparing to the region accessible for HERA kinematics one could also 
estimate that RHIC is like HERA with protons, while LHC will be more like HERA for 
(perhaps not very heavy) nuclei. 

This new range of beam energies reached at HERA, RHIC and the LHC is so high 
that "bulk" phenomena might become accessible to weak coupling calculations for two 
reasons. Firstly, given a large enough energy density in a large enough volume (in a 
nuclear collision rather than, say, a pp experiment) the sufficiently hard scale might be 
given by the temperature of the system, viewed as a blob of quark gluon plasma. Another, 
less universally accepted, conjecture is that at high enough energies, or equivalently small 
enough x, a sufficiently large momentum scale could be generated by the high density 
of virtual "wee" partons in the wavefunction of the accelerated hadron or nucleus and 
the nonlinear interactions of these partons. This latter phenomenon is referred to as 
saturation. 

The common feature in these two concepts is that they provide a possibility to 
analyze "bulk" phenomena in terms of a weak coupling constant. One could say that 
RHIC has opened up a new "firm" region 2 in the phenomenology of strong interactions, 
between the soft (hadronic and stringy) and the hard (perturbative QCD) regions. It 
could be characterized by the description of the system in terms of deconfined degrees 
of freedom, quarks and gluons, but also in terms of such high phase space densities of 
these partons that the nonlinearities of QCD must be treated nonperturbatively. As a 
digression one can mention another interesting idea prompted by the baryon excess at 
transverse momenta pt ~ 3 — 5 GeV observed at RHIC (§1 1101 II lj . This excess over 
both the thermal spectra at lower momenta and the one observed in dilute systems (pp 
collions) at the same momenta has been interpreted in terms of quark recombination 
|12j . which also relies on the large phase space density of partons. 

Because saturation is inherently a high density, or strong field, phenomenon it cannot 
be fully understood in terms of a perturbative calculation (which is a power series in the 
field strength as well as the coupling constant). However, if the saturation scale Q s 
(the momentum scale characterizing the density of partons at which the nonlinearities 
dominate) is large enough, one might still be able to perform a weak coupling calculation. 
The argument for a classical field approximation arises from these circumstances; one 
can perform nonperturbative calculations, but use the weak coupling to argue that 
quantum corrections can be neglected. The nonlinearities of the Yang-Mills Lagrangian 
(see Appendix IA.2I for the explicit form) start to dominate when both terms of the 
covariant derivative + igA^ are of equal importance. Parametrically this means that 
at saturation momentum scales, id^ ~ Q s , the gauge fields involved should be of order 
Au ~ Qs/ 9- The number density of gluons n should then be of order n ~ AA ~ Q 2 /a s . 
If the transverse phase space density is high enough, dN/ d 2 xx ~ Ql/ct s with, Q s 3> 
Aqcd; then the coupling is weak and the occupation numbers of the quantum states 
involved ~ l/a s S> 1. This is the region where one expects a classical approximation to 
be valid. 

The use of the classical field approximation in QCD is also interesting becauses it 
involves some of the most fundamental issues in modern physics; the relation between 
quantum and classical theory and the wave-particle duality. The classical field approx- 
imation is also at present the most viable practical method to study time dependent 
phenomena in field theory. 

2 I owe the term "firm" to the talk by M. Lisa at Quark Matter 2004, unpublished. 
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A classical field approximation was used to study heavy ion collisions already in e.g. 
|13| . but the idea of using the classical field approximation gained more traction with 
the model written down by McLerran and Venugopalan |14| 1151 116) . 

So far, however, most applications of saturation to phenomenology have been pertur- 
bative calculations with saturation included as some kind of phenomenological modifi- 
cation of perturbative gluon distributions E3 1211 ■ While this can be appropriate 
in some cases, such as pA-collisions at RHIC where only the nucleus is in the saturation 
regime ^3 Effl 1221 ESI or heavy quark production [24] . there are phenomena like gluon 
production in central heavy ion collisions for which one would like to have a nonper- 
turbative description. It is for this reason that also classical field theory computations 
with RHIC phenomenology in mind have been performed, starting from [25] ; they are 
the subject of this thesis. 

Now that RHIC has already been operating for several years, the best that the clas- 
sical field calculations have achieved are a posteriori explanations for what has already 
been observed. But learning from these postdictions for RHIC we should be in a posi- 
tion to make more concrete quantitative predictions for the LHC results than were made 
before RHIC operations started |2S]- For a review on applications of classical field or 
saturation ideas to RHIC phenomenology, see e.g. |2Z1 I2E1 [2Sj • 

In the following we shall first try to convey a broad picture of a relativistic heavy ion 
collision in Chapter [21 In Chapter [3] we shall discuss saturation in the proton or nucleus 
wavefunction, how it has been studied in deep inelastic scattering and how saturation 
appears in the McLerran- Venugopalan model for the wavefunction. Then in Chapter 0] 
we turn to applying these ideas to calculating particle production in heavy ion collisions. 
In Chapter [3J we discuss in more detail some of the numerical methods used in these 
calculations. In Chapter El we review the results of the three publications included in 
this thesis PJ [21 El before concluding in Chapter 
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Chapter 2 

Relativistic heavy ion 
experiments 

2.1 Spacetime picture of the collision 

We shall be interested in studying the case where two nuclei move at the speed of light 
along the = 0-axes 1 . These nuclei then collide and leave behind them, at finite 
values of i] and r, some matter which is then observed in detectors located in some 
region, varying between different experiments, around r] = 0. In order to properly 
interpret what is measured in the detectors one needs to understand the whole collision 
process. 

The common baseline for understanding the spacetime evolution of the matter formed 
in an ultrarelativistic heavy ion collision owes much to Bjorken's boost invariant hydro- 
dynamical model [3131 (see also j3J f° r a spacetime picture of the early stages of the 
collision). The Bjorken model assumes, based on experimental data from pp-collisions, 
that at high enough collision energies the central rapidity region, far enough from the 
fragmenting nuclei, can be described as a boost invariant system, i.e. with particle phase 
space densities independent of spacetime rapidity r\. This model provides a relatively 
simple framework for understanding the collision process, but contains several assump- 
tions that have to be theoretically understood and, if possible, experimentally verified, 
involving several areas physics. 

1. The initial condition at r = depends on the properties of the nuclear wavefunc- 
tion at small x and the dynamics of the partonic collision process. 

2. The thermal and chemical equilibration of the matter formed at r < To in principle 
requires understanding of time dependent, nonequilibrium Quantum Field Theory. 

3. The Quark Gluon Plasma phase lasts for some fermis during to < r < t^. For 
RHIC the hadronization timescale is < 10 fm, for the LHC it is expected to 
be larger. If the system reaches local thermal equilibrium, its behavior can be 
described using finite temperature field theory 2 . 

4. Finally, for r > 10 fm the system hadronizes and then after some time decouples 
and the particles fly to the detectors. There will be no attempt to describe in 

1 See Appendix I A . 1 1 for the coordinate system. 

2 For a review of recent developments in this field see |32|. 
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detail the setup of the detectors here, for an introduction to the detectors see e.g. 
[33] for RHIC and [MIES] for the LHC. 

The area of interest of this thesis are the first two phases of this picture, the initial 
conditions and thermalization. In Sec. l2.2l we shall discuss the more detailed experimen- 
tal evidence for thermalization in heavy ion collisions and mention more fundamental 
studies aimed at understanding thermalization as an aspect of time dependent quan- 
tum field theory. Then in Sec. 12.31 we shall turn to what can be said about these early 
stages of the collision by looking at "bulk" observables, i.e. the transverse energy and 
multiplicity. 

Apart from ion-ion collisions the nuclear wavefunction can also be independently 
studied in collisions between the nucleus and a more compact probe. The simplest such 
probe is a lepton, leading one to study deep inelastic scattering. We will discuss a 
saturation model approach to DIS in Sec. 13.21 Another possibility is to use a proton 
or, as is done at RHIC for technical reasons, a deuteron. Measurements of high px 
particle spectra and jet-like correlations at RHIC j3j3 OH EHl EH] have been essential 
in separating initial state nuclear or saturation effects, visible already in dAu-collisions, 
from effects caused by the presence of a strongly interacting medium, i.e. final state 
effects 3 . 

2.2 Thermalization 

The matter formed in the early stages of the collision is evidently so dense and strongly 
interacting that it is opaque to high pr jets [15]. This can be inferred from the suppres- 
sion of high pt hadrons in central AuAu-collisions |44| 13 7| I36j and the disappearence of 
jet-like correlations between particles emitted in opposite azimuthal directions |36l I45j . 
But is this matter in thermal equilibrium? 

The yields of different species of hadrons finally measured are remarkably well repro- 
duced by a thermal fit [M] depending only on the decoupling temperature Td ec , volume 
and baryochemical potential. This is an indication that the hadrons are emitted from 
an equilibrated system at T = Td ec . It has been noticed, however, that the same kind of 
fit also works for pp and even ee-collisions |47| 148 \ I49j and might thus result from from 
a purely combinatoric argument concerning the hadronization process |5Uj . 

A perhaps more compelling argument for this scenario of thermal equilibrium and 
hydrodynamical evolution is the success of hydro dynamical models |54| in explaining a 
variety of "momentum space" experimental observables, such as elliptic flow [HE] (see 
Fig. 12.1(1 and particle spectra |561 I57j (see Fig. I2.2j) . Hydrodynamical models have 
been somewhat less successful in explaining the "position space" sizes of the system as 
measured by HBT interferometry [S3] (this is referred to as the "HBT puzzle"). 

One must, of course, also be suspicious as to how much these successes of hydrody- 
namical is actually due to a thermal nature of the system. Although one is easily led 
to conclude from the experimental observations that the system indeed does thermalize 
early enough to justify the use of hydrodynamics, it is not yet very well understood 
theoretically what exacly is the timescale of thermalization and how it could happen 
fast enough to justify the Bjorken scenario. 

3 For a theoretical calculation interpreting these measurements in terms of the "color glass condensate" 
see e.g. Ref. gOI EH 133 
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Figure 2.1: Elliptic flow t>2 of different particle species as a function of trensverse mo- 
mentum. Figure and experimental data from Ref. 51 , hydro calculations from Ref. |52|, 




Figure 2.2: Right: Transverse spectra of identified hadrons compared to hydrody- 
namical calculations, with tq — 0.6 fm. The different theoretical curves are for different 
decoupling temperatures, Td cc = 100 GeV for the thick lines and Td ec = 165 GeV thin 
lines. The dashed lines have a different initial velocity profile. Figure from Ref. |53j . 

There is a wide literature on thermalization in numerical calculations of real time 
quantum field theory |591 I6U1 1611 1621 1631 1641 165j (see [SH] for a recent review) . These 
calculations, however, have so far been mostly limited to scalar field theory (although 
very recently also gauge theory has been studied |HZ|) and to 1+1 dimensions so that 
they are not applicable to RHIC physics. One must also point out that the geometry of 
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the initial stages of a heavy ion collision is not a static 3 dimensional space, but rather 
one that is expanding in the longitudinal direction. Indeed, the phenomenologically most 
important aspect of thermalization is not necessarily that of the thermal (exponential) 
distribution of particles in different momentum modes but the isotropization of the 
momentum distribution between the transverse and the longitudinal degrees of freedom 

[nHiinni- 

Most perturbative estimates, e.g. the bottom-up scenario [HE|, generically produce 
quite a large thermalization time, tq > 3 fm. Also approaches to chemical equili- 
bration based on kinetic rate equations do not seem to reach chemical equilibrium 
fast enough to explain RHIC particle yields [70:. On the other hand, one could ar- 
gue that if the behaviour of the system is characterized by some quite large momen- 
tum scale, i.e. the saturation scale Q s ~ 1 ... 2 GeV, thermalization could occur 
already at times tq ~ 1/Q S ~ 0.2 fm. It has been pointed out recently (see e.g. 
Ref. !S3[m[I2[I3[m[ZSllZSllZI!) that plasma instabilities could provide the rapid 
thermalization that hydro dynamical models require. 

A related question is the viscosity of a strongly coupled deconfined medium. In the 
weak coupling limit the shear viscosity is proportional to l/g A 4 . In ideal hydrodynamics 
the viscosity is neglected, so it would be important to know what the viscosity is in the 
strong coupling limit. It has been proposed that the viscosity of N = 4 super- Yang-Mills 
theory in the strong coupling limit that can be calculated the AdS/CFT correspondence 
could give some insight into this question |80j . 

In this context classical field models of the nuclear wavefunction and particle pro- 
duction can provide some insight into understanding the collision process. Although the 
actual equilibrium state of a classical classical field theory is not the correct one 5 , one 
could hope that the classical theory would give a reasonable phenomenological insight 
into the timescale of thermalization. 

2.3 Transverse energy and multiplicity 

The transverse energy and multiplicity at midrapidity are perhaps the most simple 
RHIC observables and consequently the first ones measured. They are also an example 
of the kind of quantity that cannot be calculated in traditional collinear perturbative 
QCD. Phenomenological models traditionally used to estimate this kind of quantities 
have been dominated by nonperturbative ("stringy") physics (see e.g. |HB)- The idea 
of gluon saturation at small x opens up the fascinating possibility that such quantities 
could be understood with a weak coupling calculation based on the QCD Lagrangian 
with the nonperturbative aspects of the calculation factorized into properties of the 
nuclear wavefunction. 

The difference between the multiplicities and energies at yfs = 130 GeV jH21 ESI EH 
155] and y/s = 200 GeV jSHl EH IHFJ EH! is not very large at the level of accuracy of the 
present discussion, but let us for concreteness quote the values from Ref. [EH] for the 

4 For recent weak coupling calculations see e.g. |78II79| . 

5 The classical theory exhibits is the famous Ray leigh- Jeans divergence that was one of the motivations 
for Planck's quantized theory of electromagnetic radiation. The classical equilibrium state is one where 
the energy is evenly distributed between all degrees of freedom, whereas in a quantum field theory only 
modes with energy < T are populated. 
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transverse energy in central collisions at y/s = 200 GeV 



dE T 

» 620 GeV (2.1) 

dr] 

and the ratio of the transverse energy to the charged multiplicity 

— ps 0.86 GeV. (2.2) 

The particles produced in the central rapidity region are predominantly pions (to a first 
approximation equal numbers of 7r ± and 7T°), so we can approximate that the charged 
multiplicity iV c h is 2/3 of the total multiplicity, giving 

^± « 1100. (2.3) 
dr/ 

What, then, can we infer on the properties of the system at early times from these 
numbers? 

In the Bjorken scenario the system undergoes an isentropic boost invariant expansion 
staying in thermal equilibrium from the very early time until decoupling, i.e. the entropy 
in a unit of spacetime rapidity stays constant. The entropy is directly proportional to 
the number density 6 , and thus we can directly relate the measured multiplicity to the 
initial state: 

t°L ~ i°L_ ~ noo. (2.4 

dr] dr) 

In the ideal hydro dynamical Bjorken expansion the energy per unit rapidity decreases 
with the proper time as ~ r -1 / 3 . Because the volume of the unit of spacetime ra- 
pidity is r times the transverse area, this means that the 3-dimensional energy density e 
decreases like r -4 / 3 . The energy decreases because an expanding system with a pressure 
does p dV-work, i.e. the missing energy disappears down the beampipe to larger rapidi- 
ties. It is because of this phenomenon that the most important aspect of thermalization 
for the energy and multiplicity is the creation of a longitudinal pressure |9Uj . The es- 
timate of Ref. [H5> assuming an early time for starting the decrease in the transverse 
energy, is that the it could decrease even by a factor of 3.5 between r = 0.2 fm and 
decoupling. As the assumed thermalization time is very early, this could be considered 
an upper limit, giving an estimate for the initial state: 

dE T (r = 0.2 fm) ^ . . 

- < 2200 GeV. (2.5 

dr] 

To express the energy density in units of GeV/ fm 3 that are often used one must specify 
the proper time. Let us somewhat arbitrarily consider r = 1 fm as a typical early time 
in the collision process when hydrodynamical evolution could be taken to begin. The 
decrease of a factor of 3.5 in Ref. [HJ assumed a very early thermalization time. If the 
hydrodynamical evolution is started only later at r = 1 fm, the transverse energy will 



6 The coefficient is different for bosons and fermions and depends on the masses of the particles, but 
as the measured particles are mostly relativistic pions (bosons with (jar) ~ 0.5 GeV > m n ) and the 
particles in the initial state mostly gluons, we shall neglect these factors in this discussion. 
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only be able to decrease by a factor of 3.5 x (0.2 fm/1 fm) 1//3 ~ 2. This leads to the 
estimate 

e(r = 1 fin) = ^ = 1 fm) * < 9 GeV/ fm 3 . (2.6) 
dn tttR\ 



for the 3-dimensional energy density at r = 1 fm. Hydro dynamical calculations with 
initial conditions from different saturation models have been performed in e.g. |57| I91| 

Another limit may be obtained by proceeding as in e.g. |171 1181 1201 1211 1411 I95j 
where, analoguously to e.g. pp collisions one tries to directly relate the multiplicity of 
partons in the initial state to the final multiplicity without assuming thermalization 
or work done by a transverse pressure, simply assuming that the ratio of the parton 
multiplicity in the initial state and the hadron multipliciy in the final state is constant. 
This idea is in a sense related to the so called parton hadron duality [HE] ■ If the system 
is not in thermal equilibrium, the entropy and thus the particle number density is not 
conserved but increases . It has also been pointed out that at least in a perturbative 
calculation particle number conserving collisions thermalize the partonic system quite 
slowly and that particle number increasing processes are essential for thermalization 
|68| . causing the multiplicity to increase. 

If there is very little longitudinal pdV work done by the pressure the energy in a 
unit of rapidity will be conserved: 

dE init 

—Zj- > 600 GeV, (2.7) 
meaning that the energy density at r = 1 fm would be 

J pinit -i 

e(r = i fm ) = —L > 4 GeV/ fm 3 . (2.8) 

dn tttR a 

The multiplicity, on the other hand, can increase, so our estimate for the initial multi- 
plicity, Eq. (|2.4|) should be taken as an upper limit: 

A /v init A AT final 



dn ~ dn 

Essentially the difference between these two scenarios boils down to the difference 
between the energy density decreasing as 1/r for a free streaming scenario and 1/r 4 / 3 
for isentropic ideal hydrodynamical expansion. If the system equilibrates early and 
decouples late, the difference between 1/r vs. 1/r 4 / 3 will be large. 

We will return to these estimates and their meaning in terms of the classical field 
model in Chapter El 



7 The same happens if the matter is not an ideal fluid but has viscosity. 
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Chapter 3 



Parton saturation in the small x 
wavef unction 



3.1 Parton saturation 

The gluon distribution in a hadron or a nucleus, measured at a fixed virtuality Q 2 , is 
seen in deep inelastic scattering experiments to grow towards smaller x |98[I99|. This can 
be understood as resulting from a cascade of gluons radiated from the larger x degrees 
of freedom into the larger phase space made available by the increasing energy yfs 1 . 

The general idea of parton saturation jlUUl HUlj is that as the transverse phase 
space density of gluons grows large enough, it starts to be limited by gluon fusion or, 
equivalently, the nonlinear interaction of the color field with itself. More specifically, the 
nonabelian field strength tensor consists of the linear part dA and the nonlinear part gA 2 , 
and for d ~ pt < 9^ the nonlinearities start to dominate. In this scenario the number 
density of gluons would be limited from above by n ~ AA < l/a s . The momentum 
scale below which this happens is called the saturation scale Q B (x). Different authors 
use different conventions for the saturation scale (we shall mention some in Sec. 13.4(1 . 
but to illustrate the idea let us cite one definition from |102| : 



This is formally an implicit equation for Q s (x) but because the r.h.s. depends on 
Q s (x) quite weakly (logarithmically) it is easy to get a reasonable approximation for 
the solution. In Eq. (|3.1|) p is the nuclear (baryon number) density, and together with 



2^J R\ — b|,, the longitudinal extent of the nucleus at impact parameter b^, it gives a 
momentum scale that is inversely proportional to the transverse area of the nucleus. If 
one argues by the uncertainty principle that a gluon with transverse momentum Q s has 
a transverse area ~ l/Q 2 , Eq. dSHJ) just gives an explicit definition for the saturation 
scale as the scale where there are on the average ~ l/a s gluons overlapping in the nu- 
clear wavefunction. When the strength of the interaction between these ~ l/a s gluons 
is ~ a s , this means that at the saturation scale (and lower momentum scales) the gluons 
cannot be treated as independent. 

Or actually v, the photon energy in the target rest frame, in DIS experiments. 




(3.1) 
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These saturation ideas have been implemented in various ways. One straighforward 
application in the framework of the traditional collinear factorized formulation is to ap- 
proach saturation from higher momentum scales, and include the first nonlinear correc- 
tion (the GLRMQ terms named after the authors of in the DGLAP evolution 
equations. This approach has been applied e.g. to charm production |1L)31 I1U4"] . which 
is an example of a quantity that is most naturally computed in the collinear factorized 
framework. The saturation argument has also been used to construct a successful model 
(the Golec-Biernat & Wusthoff model |1()5| 1106] ) for describing the HERA data for deep 
inelastic scattering at small x. We shall discuss this model in Sec. 13.21 

For decreasing x which, for a fixed Q 2 , means increasing energy, Q s is expected to 
grow. For large enough energy (or for large enough nuclei) the limit Q s 3> Aqcd should 
be reached, making it possible to use weak coupling methods. Because the color fields 
are strong, the occupation numbers of quantum states of the system are large, and it 
is natural to use a classical field approximation. We shall discuss the classical field 
approximation in Sec. 13.31 before briefly reviewing some other points of view on parton 
saturation in Sec. 13.41 



3.2 Saturation in DIS 




Figure 3.1: In the dipole frame the incoming virtual photon splits into a quark- 
antiquark dipole of trensverse size r^, which then interacts with the target with the 
dipole cross section a. 

It is useful to think of deep inelastic scattering at small x in the dipole picture 
|1D7[ I108[ I109| I110| llll[ I112j , where the process is viewed as a virtual quark fluctuating 
into a color dipole, which then probes the wavefunction of the target (see Fig. 13.1(1 . 

In the dipole model one factorizes the total cross section into the probability for the 
virtual photon to fluctuate into a qq pair (a color dipole) and the cross section of the 
dipole scattering with the target 2 . The total cross section can be written as |105[I112*] : 

a T ,L(x,Q 2 ) = J d 2 r T ^ dz\^ L (z,r T )\ 2 a(x,Q 2 ,r T ). (3.2) 

2 There are some tricky issues related to the Lorentz frame in which one should view the scattering 
process in the dipole model, see e.g. the discussion in |113l I114| . 
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Here the photon wave function ^t,l{ z i v t) gives the probability for the virtual photon 
(T and L stand for, respectively, transverse and longitudinal polarizations of the photon) 
to split into a color dipole of transverse size r^. The wave function ^t,l{z,yt) includes 
the known QED part of the reaction and is known analytically 3 . The exact expressions 
can be found in e.g. Ref. |1U5| . 

It has been shown PUSES] that the HERA data on the total [TTTfllTT7niTT71ITT7aiTT3] 
and the diffractive |1201 11211 1122j cross sections for x < 0.01 is well reproduced by a 
saturation parametrization 

a(x, Q\ r r ) = <r (l - e^Q-OW), (3 . 3) 

with the saturation scale Q 2 depending on x by 

Ql{x) = (x/x )- X GeV 2 . (3.4) 

The values found in [TUSl were = 23.03 mb, x = 3.04 • 10~ 4 and A = 0.288 4 . 

These fits are an example of a more general hypothesis that the behavior of the 
system is controlled by a universal saturation scale Q 2 (x). A simple demonstration of 
this idea is studied in e.g. [HI I123j . where a large variety of small x deep inelastic 
scattering data from different experiments, with both protons and nuclei |124l il25, 126 , 
is found to follow a universal curve when plotted as a function of the scaling variable 
Q 2 /Qs ~ Q 2 x x , instead of being functions of Q 2 and x separately. 

3.3 The classical color field of a high energy nucleus 

The classical equations of motion of a nonabelian gauge field theory |127j are an interest- 
ing subject in themselves. But they are useful especially because in some circumstances 
the classical field approximation of the true quantum theory can be a very useful tool 
in understanding phenomena where high particle densities are important. One example 
of such a system are the soft bosonic modes of finite temperature field theory [HSj . The 
example that concerns us in this work is the small x wavefunction of a proton or a nu- 
cleus, where, as we argued in Sec. 13.11 the density of gluons grows large (parametrically 
l/a s ) and the high density gluonic system that is born when these quasireal gluons are 
freed in a relativistic heavy ion collision. 

The McLerran-Venugopalan model |14| 1151 116j is a classical field model for the small 
x wavefunction of a hadron or a nucleus. The correspondence between the classical field 
model and a diagrammatic calculation is explored in 1281 1129] . Supplemented with the 
formulation of a collision of two ions in [1301 1131 j the McLerran-Venugopalan model 
forms the basis for what is referred to in this work as the classical field model for heavy 
ion collisions. We shall now concentrate on the wavefunction of one hadron or nucleus 
and return to colliding two of them in Chapter HJ 

The central idea behind the McLerran-Venugopalan model and one that remains at 
the heart of the more general concept of the Color Glass Condensate is the separation 
between hard (large x) and small x degrees of freedom. The former are treated as 
classical sources of radiation and the latter as a classical color field generated by these 

3 To leading order in a om , which is quite sufficient in this context. 

4 These are the results with charm quarks included. Without charm the values are somewhat different 
but arguably the most important parameter in this context, A, only changes to A = 0.277. 
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sources. The idea is thus to have a classical effective theory for the small x degrees of 
freedom. 

Let us consider a hadron or a nucleus moving at light velocity in the positive z 
direction. The large x degrees of freedom (the valence quarks, in a first approximation) 
are considered classical current 

J» = 5 fl+ 5(x~)p(x T ). (3.5) 

Let us try to justify this form somewhat. The current only has a +-component, being 
caused by particles with a large +-component momentum. The naive explanation for 
the delta function 5(x~) is that when the nucleus is traveling at the speed of light, it is 
Lorentz-contracted to an infinitesimal thickness. The more proper justification is based 
on the Heisenberg uncertainty principle. The large x degrees of freedom have a large 
p + , and are therefore well localized in x~ . The small x degrees of freedom, the gluons 
that form the classical field, have a smaller p + and are spread on a larger distance in x~ , 
effectively seeing the sources as a delta function in x~ . The source is taken to be static 
in the sense that it does not depend on the light cone time x + , because due to Lorentz 
time dilation the evolution in x + of the source is much slower than the timescales of the 
small x degrees of freedom that we want to probe. 

The classical color fields representing the small x degrees of freedom are then com- 
puted using the Yang-Mills equations of motion 

[D^F^] = r. (3.6) 

To consistently serve as the source term of the classical equations of motion the current 
(|3.5|) must be covariantly conserved 

[^,J"]=0. (3.7) 

In writing down Eq. Q3.5|) we have implicitly assumed a gauge where A~ = 0. To make 
the source Eq. (|3.5[) fully gauge invariant one would have to insert a Wilson line along 
the x + -axis 132, 133 and write the source as 



J» = 5^5{x-)W\k t , x + )p(x T )W(x T , x + ), (3.8) 
where the temporal 5 Wilson line is defined as 



W(x T ,x + ) = Pexp^ig J dy + A (x T ,y + )j 



(3.9) 



What yet remains unknown is the transverse color charge density p(xj-) in the current 
(|3.5|) . The argument of the original McLerran-Venugopalan model |14| I15| I16j , put on 
a firmer group theoretical footing in |134j . is the following. Let us assume that the 
charge density is a sum of independent color charges of a large number of hard partons. 
Then the resulting charge density at different points of the transverse plane should be 
uncorrelated. The source as a whole should be color neutral, and thus the expectation 
value of the charge (p(xr)} should be zero. One can also argue by the central limit 

5 Temporal in the sense that x + is the light cone time. 
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theorem that the distribution caused by a large number of independent charges should 
be Gaussian. One is thus lead to a Gaussian probablity distribution of charges 



(p a ( XT )p b (y T )) = 5 V(5 2 (x T -yT) or (3.10) 
(p a (k T )p b (p T )) = (2ir) 2 g 2 p 2 5 2 (k T + PT ) (3.11) 

One could also argue |135[ I136j that the charge distribution should impose color neu- 
trality at some confinement scale pt < Aqcd ■ I n momentum space this means replacing 
the correlator (|3.11j) by 

(p a (k T )p b ( PT )) = (2vr) VV 5 2 (k r + p T )/(k T ), (3.12) 

where /(k T ) « for |k T | < Aqcd and /(k T ) m 1 for |k T | > Aqcd 6 - 

One can now solve the equations of motion (|3,fi[) . The solution is most easily found 
in the covariant gauge d^Acov = 0. One can find a solution with only one component 
of the gauge field is nonzero, namely A+ ov (xt,x~). In this case Eq. (|3.6|) becomes a 
2-dimensional Poisson equation 

-V 2 T A+ y = <y(x-)p(x r ). (3.13) 

We can formally write the solution as 

At ov = -5(x~)p(x T )/V 2 T . (3.14) 

Note that there is an infrared singularity in Eq. 1)3.14(1 . The most natural prescription 
to solve this amibiguity is to impose the constraint J d 2 XT/o(xT) = 0, i.e. to require 
that the source as a whole is color neutral. Imposing color neutrality at a shorter length 
scale will also remove this ambiguity. 

The covariant gauge solution has the advantage of being localized on the light cone 
in the t, z-plane, but its interpretation in terms of partons is not very clear. To interpret 
the classical field in terms of quasi-real Weizsacker- Williams gluons we must transform 
the field into the light cone (LC) gauge. This gauge transformation can be done using 
the path ordered exponential 

U(x T ,x~) =Pexpi ig J dy~A+ ov (x T , y~) I , (3.15) 



giving 



A± c = (3.16) 
Al c = ^(x T , 2 ;-)a j [/ t (x T ,x-). (3.17) 

The light cone gauge solution is not localized on the x + -axis, unlike the one in covariant 
gauge. Instead, for x~ > it is a transverse pure gauge field. The field strength tensor 
F^ u , however, is nonzero only on the light cone x~ = 0. 

The problem of understanging the small x wavefunction of the hadron or nucleus 
has now been reduced to a simple model depending on one phenomenological parameter 



6 In numerical simulations this kind of a modification is essential to control the Coulombic "tails" of 
the classical fields when studying finite size nuclei |[TI I13TI Tl38| . 
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describing charge density of the hard sources, fj,, the gauge coupling g, which, in the 
classical field approximation, is just a constant, and the geometrical properties of the 
hadron or nucleus. One can try to estimate the value of \i from the parton distribution 
functions |139| or just treat it as a parameter to be determined from experiment. 

The source charge density [i increases as one probes smaller values of x, because the 
number of partons with higher x than the value of interest are counted in the classical 
source. To study this one needs a more refined description of the longitudinal structure of 
the source |14()| . One can then develop a renormalization group equation to see how the 
hard source develops as gluons of smaller and smaller x are integrated out of the classical 
fields and included in the source. This renormalization group equation is the so called 
JIMWLK equation [T321 ITIT1 ITiSl ITiSl ITU [TIKI IHH1 ITT71 ITlkl ITimirnTfllTimiTT^ITB!] . 



Relation to the dipole model In the classical field approximation the dipole cross 
section can be expressed in terms of Wilson lines in the background field of a nucleus 
[Tnl 11551 ITS!)! ITBTj 

a(r T ) = -^J d 2 b T Tr(|l-?7 t (b T + r r /2)[/(b T -r T /2)^) (3.18) 

where U(xt) in the Wilson line in the fundamental representation: 

U{xt) = P expUgj d X -A+ v {x T , x~)\ . (3.19) 

The correlators of the Wilson lines in the McLerran-Venugopalan model were studied 
numerically in |158| and analytically in |159| . For a discussion of the different represen- 
tations and conventions on the saturation scale see also [160 . Ref. |161j discusses the 
relation between different approaches to JIMWLK equation and how the same equa- 
tion arises when considering the quantum evolution as a property of either the target 
wavefunction or the (dipole) probe. 

By calculating the gluon distribution in the McLerran-Venugopalan model one can 
relate the strength of the color source to the saturation scale Q s defined by A. Mueller 
and Yu. Kovchegov (see e.g. Ref. |138[ IT62_) 



One must also note that there is a dependence, although only a logarithmic one, on 
an infrared cutoff in Eq. (|3.20() . In analytical calculations one often argues that the 
relevant cutoff is Aqcd> because at that scale confinement physics sets in. In numerical 
calculations the implicit infrared cutoff is the system size ttR\. 



3.4 Other views of saturation 



Let us then briefly mention some other views and aspects of parton saturation before 
proceeding to calculate gluon production in the classical field model in Chapter |1J 

An interesting proposal, discussed in e.g. in Ref. (1631 11641 11651 ITSBI 11(57] . is to relate 
parton saturation to percolation. As we have noted one can argue by the Heisenberg 
uncertainty principle that an individual parton with transverse momentum px occupies 
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an area ~ in the transverse plane. Saturation then corresponds to the situation 

where the wavefunctions of these partons begin to overlap and, because of their strong 
interactions, behave collectively. With this picture in mind one can consider saturation 
as a percolation phase transition. One would then expect to see clear signals of critical 
behavior in the system, say as a function of centrality in a heavy ion collision. 

The EKRT-model [0^ is a final state saturation model, where the saturation scale 
p sa t is determined from a geometrical saturation condition for the gluons freed in a heavy 
ion collision. The same scale p sa t serves as an infrared cutoff for the pQCD calculation 
giving the number of produced gluons in terms of usual parton distribution functions. 
Although the value of p S at is close to Q s , the saturation scale in the nuclear wavefunction, 
it is not possible to give an explicit formula relating the two due to the different concepts 
of initial and final state saturation. 

The saturation scale, or radius, Q s = 1/R S in the work of Golec-Biernat and Wusthoff 
|lf)5j or Rummukainen and Weigert |151j is the same as the saturation scale Q s of 
Kovchegov et. al. except with Ca replaced by Cf- The technical reason for this is that 
they consider correlators of Wilson lines in the fundamental representation, whereas the 
gluon distribution that Mueller and Kovchegov consider involves the same correlator in 
the adjoint representation. For a comment on the relation to the saturation scale used 
by E. Iancu et. al., see Ref. 

We are using the old notation of of Krasnitz et. al. |25| . In their newer work |138j 
they use the parameter A s defined as A s = g 2 ji. 
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Chapter 4 



Particle production in the 
classical field model 

4.1 The classical field model 




Figure 4.1: Spacetime structure of the KMcLW model. The gauge field in regions (1) 
and (2) are pure gauge fields on one nucleus. They are used to find the initial condition 
to find the gauge field in the region (3). 

Let us then turn, following the approach of Kovner, McLerran and Weigert (KMcLW) 
|131) . to studying the collision between two nuclei in the McLerran- Venugopalan model. 
We want to calculate the classical gauge field in the forward light cone using as initial 
condition the fields of the two nuclei, given by Eq. (|H.lti|) . This calculation in the light 
cone gauge was formulated and performed perturbatively to leading nontrivial order in 
|LS()[ ILS1[ ll.'fflj and to leading order on one of the sources and all orders in the other in 
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|168| . The same computation can also be done in the covariant gauge |169[ll7fH 1162 j. 
We start with a static current of the two nuclei, 

J" = ^ + P(i)(x r )<5(x-) +^-p (2) (x r )5(x+), (4.1) 

anticipating the fact that we will be working in a gauge where this form is covariantly 
conserved. In the regions x~ > 0, x + < (1) and x + > 0, x~ < (2) (see Fig. 14.1(1 the 
field is given by the pure gauges as in Eq. (|3.16j) 

A U = l -e iA ^d ie - iA ^, with V^A (m) (x T ) = - 5 p (m) (x r ), m = 1,2. (4.2) 

Here we have introduced the notation A( m ) for the solution of the Poisson equation. 
The quantity A( m ) is related to to the covariant gauge fields of the nuclei by gA+ ov = 
5(x-)A {1) , gA~ v = S(x+)A {2) , 

We then choose to work in a temporal gauge A T = (x + A~ +x~A + )/t = 1 . This will 
enable us to use a Hamiltonian formalism. It also matches smoothly to the conditions 
A~ = for = and A + = for x + = that, as discussed in Sec. 13, 3| enable us 
to drop the temporal Wilson lines from the current, Eq. (|3.8|) and use the simple form 
(|4.1|) . In this gauge the remaining components of the gauge field are the transverse 
components A i and the "longitudinal" component A v = —t 2 A v = x + A~~ — x~A + . 

Inside the future light cone (region (3) in Fig. 14. lj) the gauge fields satisfy the equa- 
tions of motion in vacuum. What is needed is the initial condition for solving these 
equations. These initial conditions can be obtained by requiring that the fields in the 
different regions match smoothly on the light cone. In practice the initial conditions can 
be obtained by inserting the ansatz 

A i = 9(-x + )9(x-)A l {1) + 6(x + )8(-x-)A\ 2) + 8(x + )6(x-)Ai 3) (4.3) 
A v = e(x + )6(x-)Aj 3) (4.4) 

into the equation of motion (|3.6j) and requiring that the singular terms arising from the 
derivatives of the ^-functions cancel. In this way one gets the following initial conditions 
for the gauge field in the future light cone: 

A{ 3) | T=0 = 4i)+4o ( 45 ) 
4^=0 = |[4),4)]. (4.6) 

The equations of motion with these initial conditions can then be solved either numer- 
ically or perturbatively in the weak field limit. Let us first look at the equations of 
motion in more detail and then review the weak field solution of |131j in Sec. 14.31 

4.2 Boost invariant Hamiltonian 

The initial conditions, Eqs. (|4.5|) and (|4.6j) . are boost invariant. Therefore it is natural to 
assume that also the solution of the equations of motion will be independent of rapidity. 
To maintain a consistent boost invariance of the solution one must not perform gauge 
transformations that depend on the rapidity rj. This limitation reduces the longitudinal 

1 See Appendix IA. II for our conventions concerning the coordinate system. 
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gauge field to an adjoint scalar field, and we will denote it by A v = cj). We could 
equally well choose = —A^/t 2 as our canonical variable; this would simply cause 
powers of r to appear in different locations (see Appendix IA.1|) . The appearence of an 
adjoint scalar is analoguous to the way the time component of the gauge field becomes 
an adjoint scalar when dimensionally reducing high temperature gauge theory to a 3 
dimensional effective theory jo2j . Gauge field theory with an adjoint scalar field in the 
r, ^-coordinate system has been extensively studied in [T7T1 1X72111751 IT7I1IT751 IT7S1 IT77| . 

In the gauge A T = and with the assumption of boost invariance we end up with 
the following action of a 2+1-dimensional gauge theory with a scalar field in the adjoint 
representation, 

S = J dr/d 2 x r drrTr|i4i-^i^ + ^ 2 -^[A^][A,^]|. (4.7) 

The fields in this action are now all functions of the proper time r and the transverse 
coordinate xy. A dot denotes a derivative with respect to r, i.e. Ai = d T Ai. The 
explicit time dependence in the action is caused by the r, ^-coordinate system. The 
expanding longitudinal geometry of the coordinate system reflects the physical situation; 
the collision region is expanding with light velocity as the two sources move apart from 
each other. 

In the A T = gauge we can go to the Hamiltonian formulation easily. Defining the 
canonical momenta (electric fields) 

E ia _ = T ia = _ TjL ai ( 4 . 8 ) 

TT a = ^- = -4> a (4.9) 
5(p a r 

we get the Hamiltonian density |172j 



H = 2Tr[£Mi + 7T0]-£ (4.10) 
Tr + ^FijFij + rvr 2 + i[A, 0][A 

The first term in the Hamiltonian is the kinetic energy term of the transverse electric 
fields. The second term is the potential energy of the transverse gauge fields, which, in 
a 3-dimensional language, would be called the z-component of the magnetic field. In 
2+1-dimensional gauge theory the magnetic field only has this one component. The 
third and fourth terms are the kinetic and potential energy, or electric and magnetic 
terms, of the scalar field. 

The Hamiltonian equations of motion in the vacuum can be obtained in the standard 
way: 

r it i 

Ai = t a j- = -E* (4.11) 
5EI t 

<f) = t a ^— =TTT (4.12) 

Y 5p a K J 

& = -t a ^=r[D kl F kl ]- l j-[^[D il 4>]] (4.13) 

* = - <a ^ = ^ D «.[A,# (4-14) 



19 



Initial conditions for the Hamiltonian variables The initial condition for the 
transverse gauge fields is given directly by the sum of the two transverse pure gauges, 
Eq. (|4.5|) . Because the transverse electric field is proportional to r: E l = rA{ its initial 
condition is simply E % [j = 0, Xt) = 0. The initial condition for A r] is given by the 
commutator of the two transverse pure gauges, Eq. Because (f) = A v = -t 2 A^, 

this means that (f>(r = 0,xy) = 0. The corresponding momentum, on the other hand, 
is 7r = 4>/t = —2A V — tA v and thus the initial condition for ir is tt(t = 0,x^) = 



-w 



4* A i 



4.3 Gluon production in the weak field limit 

Let us then calculate the number of gluons produced in a collision of two ultrarelativistic 
nuclei to the leading nontrivial order in the classical sources. Our calculation is essen- 
tially the same as performed in |131j , exept that we will use the notation of the previous 
section. We will expand in powers of the dimensionless functions A^^x^) defined in 
Eq. (j4.2|) . They are proportional to the source strengths, A( m )(xy) ~ gp^ixrr) ~ <? 2 /U, 
so the dimensionless parameter that must be small in this computation is g 2 fi/k.T- We 
will have to expand the fields to order A 2 and thus the multiplicity, which is quadratic 
in the fields, will be proportional to A 4 . 

Before solving the equations let us first discuss the definition of the multiplicity. Our 
starting point is the requirement that the energy can be expressed as an integral over 
momentum modes: 

d 2 k T n(k T )|k T |. (4.15) 



H 



This defines the differetial multiplicity n(k^), assuming a massless dispersion relation 
£"k T = |kr| and a given partition of the energy into momentum modes. We shall also take 
advantage of the equipartition of energy between the coordinates and the momenta in a 
classical mechanical system and only use only the momentum part of the Hamiltonian 
from Eq. (|4.10j) . 



H 



2 I d z x T Ti 



d 2 k 7 



(2vr) 2 



Tr 



-E i E i + T7T 2 

r 



(4.16) 



-^(kT^X-kr) + rvr(k T )vr(-k T ) 



to define the multiplicity of gluons. The "~" should be taken as an equality in a time 
averaged sense. The advantage of using only the momenta is that we can easily express 
the kinetic part of the Hamiltonian in terms of the Fourier transforms of the momenta. 
Equating the expressions ()4.15|) and (|4.16f) we then find the differential multiplicity 



n(k 



1 



T) 



(2tt) 2 |k T | 



-^(k r )S*(-k T ) +r7r(k r )7r(-k T ) 



(4.17) 



Unlike the energy, the multiplicity as defined in Eq. (|4.17|) is not a gauge invariant 
concept. The prescription we shall use is to fix the Coulomb gauge in the transverse 
plane, diAi = 0. 
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The initial conditions for the canonical variables are, expanded to order A 2 , 

Ai(0,x T ) = — (Q(A (1) (x T ) + 3;A (2) (x T )) + [<%A (1) (x T ), A (1) (x T )] 
9 ^9 

+^[d i A (2) (x r ),A (2) (x T )] (4.18) 
vr(0,x r ) = - [5 i A (1) (x r ),5 i A (2) (x r )] . (4.19) 

It is easiest to fix the Coulomb gauge already in the initial condition. This removes the 
lowest order terms in A^r = 0,xt), giving 

4 C ° ul (0,x T ) = ± fa - ||) ( [A (1) (x T ),a i A (2) (x T )] 

+ [A (a) (x T ),a 7 -A (1) (x T )]Y (4.20) 

From here on we will drop the superscript "Coul" and consider all fields in the Coulomb 
gauge for the rest of this section. Now that both A{ and tt are of order A 2 , it is easy to 
linearize the equations of motion Eqs. I|4.11|) — (|4.14|) . They can then be solved by Fourier 
transforming with respect to the transverse coordinate 

& = d T (TAi) = rV 2 T Ai (4.21) 
=> (T 2 d 2 +rd T + T 2 k T 2 )Mr,k T ) =0 (4.22) 

tt = d T (J-<fi\ = iv|0 (4.23) 

(r 2 d 2 - rd T + r 2 k T 2 ) 0(t, k r ) = 0. (4.24) 

The solutions of these equations are Bessel functions 

MT,k T ) = ^(0,k T )J (|k T |r) (4.25) 
0(r,k r ) = -1-^(0, k T )J 1 (|k T |r). (4.26) 

Using the asymptotic expansions of the Bessel functions we get the expectation value 
of the multiplicity 

(n(r,k T )> = ^^{k^sin 2 (|k T |<r- J) (A«(0, k T )A?(0, -k T )) 



(2vr) 2 vrk T 2 



+ sin 2 (|k T |r-^ ) (-"(0.k 7 )^'(0.-k 7 )) } ( 1.27) 



The correlators (A£(0, k T )Af (0, -k T )) and (vr a (0, k T )vr a (0, -k T )) can be calculated us- 
ing the initial conditions, Eqs. (|4.18f) and (|4. 19f) . the Poisson equation (|4.2|) relating the 
A's to the charge density p and the correlator of the charge densities, Eq. ()3.1()j) . One 
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obtains 



(Af(0,k T K(0,-k r )) = -R\ Nc{N y l \ 9 2 ^ f (4.28) 



,2/, xPt 2 Qt 2 - (pt • Qt) 2 

(5 (k r - pt - qr) — 2 — i — 5 

qr Pr 

2 i\ /• J2 j2. 



<vr«(0,k r K(0,-k T )> = v tf A M^ tfrf J ±P0^ (4.29) 

(pt • qr) 2 



6 (kr - pt - qr) 



q T 4 P T 4 



Because we only used the momenta to define the multiplicity, we are left with the 
oscillating sin 2 |ky|r factor in Eq. (|4.27|) . We will replace this oscillating factor with its 
time average of 1/2, which is equivalent to using also the contribution from the fields. 
To be more explicit; when we argued by the equipartition of energy and chose to use 
only the momenta, we effectively replaced 1 = sin 2 |kr|r + cos 2 |kr|r by 2 sin 2 |kr|r. If 
we now use the time average sin 2 (|kr|r) = 1/2 the result will be equivalent to using 
both the fields and the momenta. 

With this prescription the (py • qT) 2 -terms cancel between the transverse and scalar 
fields and we end up with the infrared divergent result 



(2i) 2 k T 2 ii (2i) 2 p T 2 (k T - Pt) 2 



This is the result found in |131| (some small errors of |131j were corrected in |139j , where 
the constant factors are correct). This result is analoguous to the bremsstrahlung result 
of Gunion and Bertsch |178j but with a different, in our case infrared divergent, form 
factor for the source (see discussion in |139j ). 

If we regulate the integral Eq. Q4.3UJI with a mass m, replacing p^ 2 and (ky — Pt) 2 
with pt 2 + m 2 and (kr — pt) 2 + m 2 , it can be evaluated and we obtain 



ttR\ 1 iV c (iV c 2 - l)gV kr 2 
(2vr) 3 vr k T 4 R m 2 



{n(k T )) = 7 ^-- — hi — . (4.31) 



This result still has an unintegrable singularity at ky = and therefore does not give 
a finite result for the total multiplicity. It nevertheless has some important properties 
that will also hold for the full numerically computed result. The result is proportional 
to the transverse area of one nucleus, itR\. If we assume that the infrared divergence 
if the integral (|4.3U|) should be regulated at the scale m ~ g 2 fi, where, as discussed 
previously, our perturbative expansion breaks down, the multiplicity can be written in 
a general form as 

<n(k T )> = -^f(k T /g 2 v). (4.32) 

g 2 

The numerical calculation shows that when all orders in the source are included in 
the calculation the integrated multiplicity is indeed finite. The total multiplicity is then 
given in the form 

dN J d 2 k T (n(k T )> = ^(g 2 v) 2 f N , (4.33) 



dr) 
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where 




(4.34) 



is a numerical constant (i.e. independent of g,irR\, and fj,). By the same kind of 
argument the transverse energy should be given by 

^ = J d 2 k T (n(k r ))|k T | = ^{g^ff E . (4.35) 

Numerically computing the constants fisr and Je is the purpose of the first article in 
this thesis 1 . We will discuss the numerical method for this calculation, developed in 
|25j . in Chapter |21 and turn now to quark pair production. 



4.4 Fermion pair production 

Given the classical fields corresponding to gluon production a natural question to ask 
is: how are quark antiquark pairs produced by these color fields? Formally quark 
production is suppressed by a s and group theory factors compared to gluons, so in a first 
approximation we should be able to treat the quarks as a small perturbation and neglect 
their backreaction on the color fields (see e.g. |179j for a comparison of quark antiquark 
pair and gluon production in a pQCD framework). Heavy quark production is calculable 
already in perturbation theory and in the weak field limit quark pair production from 
the classical field model reduces to a known result in /cr-factorized perturbation theory 
|180j . It is, however, unknown how much the full inclusion of the strong color fields 
changes the result. Understanding light quark production would address the question 
of chemical equilibration; turning the color glass condensate into a quark gluon plasma. 
There have been calculations of quark production in heavy ion collisions from instantons 
|181j . from static "stringy" longitudinal electric fields using the Schwinger mechanism 
[IS21 HHSl HE1 CHS1 and in the pQCD + saturation model [HD HZOJ- The approach of 
calculating the fermion Green's functions in a space-time dependent background gauge 
field in Refs. |186| I187| is more similar to our approach in Ref. |3j, but has not been 
directly applied to the gauge fields from the classical field model discussed in Sec. 14.11 

The calculation of quark pair production, outlined in more detail in Ref. [2], proceeds 
by solving the Dirac equation in the background color field of the two nuclei. This can be 
done analytically for the Abelian theory. The QED calculation, of interest for electron 
positron pair production in ultraperipheral collisions, is done e.g. in Ref. 188 . 

Let us briefly summarize the theoretical background of the calculation developed 
in Ref. |189j . It is well known that the cross section to produce one quark antiquark 
pair can be calculated by evaluating the Feynman, or time ordered, propagator for the 
spinor field. The time ordered propagator is also the quantity that is usually obtained 
when computing Feynman diagrams, because the Wick theorem applies to time ordered 
propagators. The main observation of Ref. |189j was that the expectation value of the 
number of quark antiquark pairs, is related to the retarded quark propagator. Because 
the Wick theorem does not apply to the retarded propagator it is difficult to calculate in 
perturbation theory. However, because calculating the retarded propagator is equivalent 
to solving the Dirac equation as an initial value problem, it is easy to formulate as a 
numerical calculation. One has to solve the Dirac equation with a negative energy 
plane wave as the initial condition. The resulting spinor is then projected to positive 



23 




Figure 4.2: The quark antiquark pair production calculation in spacetime. We start at 
x + < 0, x~ < with a negative energy spinor on shell, e iq ' x v(q). When the antiquark 
line meets the first color source, it is kicked off mass shell by the color field of the source 
in the form of the Wilson line J7m or U^)- When colliding with the second color source 
the antiquark gets another kick. The label U(rJJn\ refers to the Abelian case, in which 
the field in the future light cone is a pure gauge field given by the product UruUi^y In 
the non- Abelian case the field is known only numerically. In the end the wave function 
is projected to positive energy states e~ lp ' x u(p) to find the pair production amplitude. 



energy states at long enough times after the collision to obtain the quark antiquark pair 
multiplicity. 

One can give a physical interpretation for this calculation in terms of the Dirac hole 
theory, where the vacuum is interpreted as a Dirac sea with all the negative energy 
quark states filled. Pair creation by an external field is then interpreted as the gauge 
field giving enough energy to a negative energy quark to lift it to a positive energy 
state, leaving a "hole" in the Dirac sea of negative energy quarks. This hole is then 
interpreted as an antiquark. The Dirac sea interpretation should of course be considered 
as a pedagogical explanation. The real justification for this approach is the detailed 
calculation in Ref. |189j . where one starts by expressing the expectation value of the 
number of quark antiquark pairs in terms of creation and annihilation operators and 
then, by lengthy manipulations, relates this expression to the retarded propagator. 

The initial condition for t — > — oo is a negative energy plane wave. Similarly to the 
QED case, one can find analytically the solution for the regions > 0, x T < 0, because 
the gauge fields in these regions are pure gauges. These then give the initial condition 
at r = for a numerical solution of the Dirac equation in the forward light cone r > 0. 
To find the number of quark pairs one then projects the numerically calculated wave 
function to a positive energy plane wave at some sufficiently large time r. The structure 
of the calculation in different regions in spacetime is demonstrated in Fig. 14.21 
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A major technical challenge in this calculation is the coordinate system. In order 
to include the hard sources of the color fields, the colliding nuclei, only in the initial 
condition of the numerical calculation, one wants to use the proper time r instead of 
the Minkowski time t. Unlike the gluon production case, where one was able to assume 
strict boost invariance, one now has a nontrivial correlation between the rapidities of the 
quark and the antiquark. Thus, although the background gauge field is boost invariant, 
one must solve the Dirac equation 3+1 dimensions. 

A crucial point in the calculation is the choice of the longitudinal variable to be 
used with r. Possible choices are rj,z, or x^. To obtain the correct result we must 
have a dimensionful longitudinal variable, such as z or x^ to parametrize the r = 
surface. This is because one must be able to represent longitudinal momentum scales 
in coordinate space. For r > the corresponding longitudinal coordinate could be 
constructed as re _,? , but for r = this is not possible. To enable a symmetric treatment 
of both branches in Fig. 14. 2\ one should choose z as the longitudinal variable, with 
x ± = (Vt 2 + z 2 ± z)/y/h = (\z\ ± z)/V2 at r = 0. 

One is thus led to solve numerically the Dirac equation in a curvilinear coordinate 
system. The covariant formulation of spinors in a curved space is a fascinating subject in 
itself [1901, but as in this context one is merely decribing ordinary Minkowski spacetime 
with curved coordinates it was ultimately found in Ref. [2] that the most practical 
solution for this particular case is to simply change the coordinates without transforming 
the spinors themselves. The Dirac equation thus obtained is naturally linear in the 
fermion field ip, but also has an unpleasant property. Due to the choice of coordinate 
system and the decision to use the flat vierbein the coefficients of the time and space 
derivatives in the resulting Dirac equation, Eq. (|5.25[) . depend on the coordinates. (The 
vierbein determines the relation between the flat tangent space where the Dirac matrices 
are defined to the point in spacetime. See |19()j for the theory or Appendix A of Ref. [3] 
for a brief explanation.) As will be explained in Sec. 15.21 this makes finding a stable 
discretization scheme harder, and one has to use an implicit discretization scheme. 



25 



Chapter 5 



Details of the numerical 
calculation 



5.1 Classical chromodynamics on the lattice 

Let us first review, following |25| . how the 2+1-dimensional Hamiltonian developed in 
Sec. 14.21 can be cast in a form suitable for numerical simulation on a transverse lattice 
with a finite lattice spacing. This lattice version of the calculation has been used in 
many numerical studies, e.g. [U ESQ E221 EH E21 EE] ■ 

Wilson's lattice gauge field theory was first cast in the Hamiltonian form by Kogut 
and Susskind |196| I197j . The classical Hamiltonian formulation has been used ex- 
tensively in e.g. calculations of the sphaleron transition rate in electroweak theory 
|1981 11991 12UU1 12UT] . Our notations for lattice gauge theory, formulated in terms of link 
matrices and plaquettes, are explained in Appendix I A. 31 

We will start from the Wilson action (with real Minkowski time, because we want 
to obtain the equations of motion in real time) . Then we will take the continuum limit 
in the t and z-directions, leaving only the transverse coordinate discretized. This will 
allow us to change the coordinates to r, rj and restrict ourselves to r/-independent field 
configurations. Finally we will discretize the new time coordinate r using the leapfrog 
algorithm. 

The Wilson action for gauge fields is 



where the — sign between the two terms is needed because the fields are in real, and not 
imaginary time. The coordinate x stands for the different points in the 3+1-dimensional 
spacetime lattice. We then take the continuous limit in the t and z-directions to obtain 



where we have set the lattice spacing a to unity. Here we have defined the partially 
continuous limit (continuous in the /x-direction, fj, = 0, 3 but discrete in the i-direction, 



S 




(5.1) 




1 - -^Re Tr U 1>2 j + \ £ Re Tr[M 0i -M 3i ] + Tr F 2 3 , (5.2) 
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2 = 1,2) of the plaquette 



M^(x T ) = 9 - [Aj(x r ) + 4(x T + i T )] - ^(x T )(^^(x T ))C//(x T ) 

- igAifa + iT)(^C//(x r ))C/ l (x T ) - i(^[/t( Xr ))^( XT ) 

- g 2 A^{x T )Ui(x T )A^(x T + i T )C/t(x T ) (no sum over i), (5.3) 

where the t, z arguments have been left out for brevity and It denotes a unit vector in 
the z-direction, i = 1,2. 

To enable a Hamiltonian formulation we must then choose a temporal gauge con- 
dition A T = 0. We shall also assume that the field configurations are boost invariant, 
i.e. independent of r\. To consistently impose this boost invariance we must also forbid 
77-dependent gauge transformations, which reduces the longitudinal gauge field into 
a scalar field transforming in the adjoint representation of the gauge group. We shall 
consequently denote A v = <p. With these restrictions the action becomes 

S I dvdr > W — > Tr [ Ui'Ui ) - ^ I 1 - — Uh<- I rr 



/ d^dr^r/l^Tr^M 



9 2 V N c 

^^ 2 -^£^(*-&) 2 }> (5 - 4) 



where Ui denotes the derivative with respect to r. We have also introduced the parallel 
transported scalar field 

4>i(xt) = Ui(x T )4>(x T + 'it)uJ(x t ) (no sum over i). (5.5) 

The difference (j> — <j>i is a covariant difference, because the link matrices Ui have been 
used to parallel transport 4>{xt + '\t) so that it gauge transforms at the point xy, i.e. as 

<^(x T ) =>- F(x r )C/ i (x r )y t (x T + i T ) y(x T + i T )0(x T + i T )F t (x T + i r ) x 

V(x T + i T )u}(x T )V\x T )\ =V(x T )4 >i (x T )V\x T ) (5.6) 

and can thus be compared with ^(xy) in a gauge invariant way. 

To perform the Legendre transformation from the Lagrangian to the Hamiltonian 
formalism we need to replace the time derivatives in the action (|5.4j) by the canonical 
momenta, the left invariant transverse electric fields 

2ir ■ f 

E l = —^-Trt a UiUi (no sum over i) (5-7) 
T 

and the longitudinal electric fields 

7T = <j)/T. (5.8) 

In terms of these canonical variables the Hamiltonian density becomes 

H = —TrE i E i + ^-^(l-^ r ReTrU 1>2 S \ + rTr^ 2 + i J^Tr (V - ^ . (5.9) 



9 2 V N, 



c 
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The dynamics of the theory is defined by the Hamiltonian 1)5. 9 Jl and the Poisson 
brackets between the canonical variables: 



{n a ,0 b } = 5 ab (5.10) 
{E l a ,Uj} = it a 6jUj (no sum over j) (5-11) 
{El El} = 5^/^ (no sum over j). (5.12) 

The first one is the usual Poisson bracket between two canonically conjugate variables. 
The second, Eq. (|5.11j) tells us that the transverse electric field E % generates left transla- 
tions on the SU(3) group manifold of the link matrices Ui. It is actually more proper to 
consider the Poisson bracket Eq. ([5.11)1 as the definition of the electric field E l . The re- 
lation between E l and the time derivative of the link matrix, Eq. ()5.7|) then follows from 
this definition as the equation of motion for Ui. The third Poisson bracket, Eq. ()5.12|) . 
follows directly from the second using the antisymmetry of the Poisson bracket and the 
Jacobi identity: 

{{Ei,4},U k } = -{{Ei,U k },Ei}-{{U k ,Ei},Ei} (5.13) 
= S tj fabc{E 3 c ,U k } (no sum over j). 

The equations of motion for the fields can be found by taking the Poisson brackets 
between the fields and the Hamiltonian; the equation of motion for any quantity v is 
v = {7i,v}. The equations of motion thus obtained are [TT I25j 

q 2 ■ 

Ui = i—E l Ui (no sum over i = 1,2) (5.14) 
r 

(j) = T7T (5.15) 

E 1 = IL [U li2 + U!- 2 - h.c.]- trace + -[4>i,4>] (5.16) 
2g 2 t 



E 2 = ^ [U 2 ,i + U 2 ,-i - h.c. ] - trace + ' 



7T 



^[^ + ^-2^, (5.17) 



where "—trace" means that the part proportional to the identity matrix must be sub- 
tracted, because the electric fields E % are traceless matrices. It is easy to check that these 
equations are invariant under gauge transformations that depend only on the transverse 
coordinates. Note that although the r.h.s. of the electric field equation of motion, 
Eq. ()5.16j) depends on different plaquettes, these are all based on the same point, which 
makes the equation gauge invariant. 



Initial conditions on the lattice In addition to the equations of motion we also 
need the initial conditions corresponding to Eqs. (|4.5|) and ()4.6|) in the continuum. The 
initial condition for the link matrices U^ can be written as 



Tr 



to 



(u^ + un(i+un-^- 



r(2) 



t(3)v 



0, 



(5.18) 



(1 2) 

where U- ' are the pure gauge fields corresponding to the individual nuclei, given by 



f/^ m) (x T ) = e iA (m) (x T ) e -;A (m) (x T +i T ) 5 with V 2^ A 



(m) 



-9P(m)- 



(5.19) 
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By expanding the link matrices in the limit a — > it can be verified that in the con- 
tinuum limit this reduces to Eq. (|4.5jl . Equation (|5.18j) is a very nonlinear equation for 
C/(3)i) and must be solved with some kind of iterative process. When approaching the 
continuum limit the link matrices are closer to the identity matrix and solving this equa- 
tion numerically becomes simpler. If Eq. q5.18|) is gauge transformed (with the gauge 
parameter depending on the transverse coordinate), the identity matrix in 1 + U- 
becomes a pure gauge matrix. This can be interpreted 1 so that the identity matrix 
corresponds to the gauge field in the past light cone x + < 0, x~ < 0. In this sense 
Eq. (|5.18|) is not gauge invariant, it applies in the specific, physically well motivated, 
gauge where the field vanishes in the region of spacetime where neither of the colliding 
nuclei has yet passed. 

The initial condition for the longitudinal electric field can be written as 



c/r j (xT) - 1) [ur^T) - ur^ T )) - h. c . 

. (5.20) 



2g 2 

+ (c/, t{3) (x T - i T ) - l) (t/f } (x T - ir) - ^(xt - h) ) h.c 



It is easily verified that in the continuum limit Eq. 1)5.20(1 reduces to the corresponding 
commutator equation in the continuum case, Eq. (|4.6j) . Once the more complicated ini- 
tial condition for the transverse fields, Eq. (|5.18[) . has been solved, evaluating Eq. (|5.20|) 
is a straightforward calculation. 

As in the continuum case the initial conditions for the transverse electric fields and 
the adjoint scalar field (ft are E l (r = 0) = <J)(t = 0) = 0. 



The leapfrog algorithm The leapfrog algorithm is commonly used in Hamiltonian 
time evolution of classical gauge theory. It is suited for partial differential equations that 
are second order in time or, equivalently, for systems where the fields (coordinates) and 
their time derivatives (momenta) are independent variables. The leapfrog algorithm is 
time reversal invariant, which guarantees second order accuracy in time. The essential 
idea of this algorithm is that if the coordinates are defined at times r, r+ dr, r+2 dr, . . . , 
the momenta are defined at timesteps r + | dr, r + § dr, r + 1 dr, . . . The momentum 
at r + | dr is then used to "leap" the coordinate from r to r + dr, hence the name. The 
only peculiarity particular about the equations of motion in this case, Eqs. I|5.14|l — (|5.17j) . 
is the explicit time dependence, which must be treated properly in order to maintain 
the time reversal symmetry of the algorithm. This simply means that when stepping 
the fields from r to r + dr the explicit time argument must be taken to be r + | dr. 

lr The interpretation is evident from the derivation in |25|. 
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Explicitly, the timesteps can be written as 

ig 2 dr 



Uiir+dr) 
4>{t + dr) 
E 1 {t+ dr) 



E 2 (t + dr) 



exp 



E l {r+ dr/2) £/;(r) 



r + dr/2' 
r) + (r + dr/2) drvr(r + dr/2) 



^ 1(r) + i (r + dr/2)dr 



2g 2 



U\ 2 + U\ -2 — h.c. — trace 



+ 



i dr 



r + dr/2 
E 2 (r) + 
i dr 



r+dr/2 



2 i(r+ dr/2)dr 



25 2 



^2 i + U2 -1 — h.c. — trace 



+ 



r + dr/2 



7r(r + dr) = tt(t) + 



dr 



r + dr/2 ^- 



dr/2 



-dr/2 



(5.21) 
(5.22) 

T+dr/2 

(5.23) 

r+dr/2 

(5.24) 



where the parallel transported scalar field was defined in Eq. (|5.5|) and the index 
i = 1, 2 is not summed in Eq. (|5.21|) . 



5.2 Discretizing the Dirac equation in curved coordinates 

Implicit discretization Let us first discuss the general difference between an implicit 
and an explicit discretization 202 . Consider a partial differential equation of first order 
in time 2 . Denote the values of the unknown function at a timestep n by £ n , where £ n is 
a vector whose components are the values of the function at different points in space, 
labeled by, say, j. The partial differential equation can be written as £ = D£, where D 
is some differential operator, i.e. a matrix in position space Djji. This equation can be 
discretized as Cn+i = (1 + dtD)£ n . This is an explicit discretization, because only the 
spatial derivative of the known function £ n at the previous timestep is needed. If might 
happen that the matrix (1 + dtD)jji = djji + dtDjji has an eigenvalue with absolute 
value > 1. In this case the discretization is unstable, even if the original equation is not. 
This normally gives an upper limit, the Courant-Friedrichs-Lewy (CFL) condition, for 
the size of the timestep dt. An example of an implicit discretization, on the other hand, 
would be (1— dt-D)£ n+ i = £ n , requiring one to solve a system of equations (a linear one if 
D is independent of the £) to find £ n +i at each timestep. The two discretization schemes 
are formally equally accurate (to first order in dt), but the implicit one is usually stable 
for reasonable differential operators D. 

We shall now proceed to discretizing the Dirac equation in r, ^-coordinates implic- 
itly and showing in detail how the resulting linear system can be solved using LU- 
decomposition |2U2| . Let us note, in passing, that it might well turn out that to solve 
the 3+1-dimensional Yang-Mills equations in a coordinate system where the time vari- 
able is r, one will probably also have to use an implicit discretization also. This is more 
difficult than the present case, because unlike the Dirac equation, the equations are not 
linear. Nevertheless, an implicit scheme has already been used in Yang-Mills theory 
203 , so the task should not be impossible. 



2 This is quite general, since e.g. a second order equation is equivalent to a system of first order 
equations. 
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Dirac equation in r, z-coordinates It is useful to separate the Dirac spinor tp into 
eigenvectors of 7°7 3 , by the projection operators P^ = ^(1 ± 7°7 3 )- In the t,z- 
coordinate system and for these components separately the Dirac equation with an 
external gauge field is: 



d T ip 



± 



Vr 2 + z 2 ± z 



(=f9 z V ± + ^7°(^7t ' D T - m)^ T ) =F i—ip~ 



(5.25) 



Here we are using the same gauge as in the numerical calculation of gluon production, 
namely A T = 0. We also assume that the background gauge field is independent of 
rapidity and thus the longitudinal gauge field is an adjoint scalar field. From now on 



we absorb the coupling constant into the field and denote 



gA„. To shorten the 



notation we shall also denote the timestep with a subscript n and the lattice site in the 
longitudinal z-direction with a superscript j. We discretize the Dirac equation implicitly 

as 



1 



2dr 



n+l 1 "rn-l 
\Jt 2 + z 2 ± z 



= =F 
+ i 



4dzr 
Vt 2 + z 2 ± z 



7 °(nT • Dt - m)P^ n =F i^P ± [Vi+l + <-iJ ( 5 - 26 ) 



Multiplying by 2 dr and solving we get 



l±i- 



dr 



\/r 2 + z 2 ± z 



2&ZT 



dr 



P ± t 



1=F»- 



dr 



Vr 2 + z 2 ± z 
T — dr 



n - x 1 2dzr 
Vr 2 + z 2 ± z 



2dr7 U (7 T -D T + im)P^ (5.27) 



Now we must address the question of the boundaries in the z-direction. Periodic bound- 
ary conditions can be used in the transverse direction, but not for the longitudinal 
coordinate, because the equation is not invariant under translations in the z-direction 
at fixed proper time r. We choose free boundary conditions, not setting any restriction 
on tp at the boundary. In practice this means the following. Let us denote the endpoints 
of the lattice by j = ± J, i.e. the longitudinal lattice consists of 2 J + 1 points. When 
for the points in the interior of the lattice the longitudinal derivative was discretized 
using the difference tpi +1 — tp 3 ~ 2dz(d z tpy , for the endpoints this must be replaced 
by 3ip J + ip J ~ 2 - 4V> J_1 ~ 2&z(d z ip) J and 4ip~ J+1 - 3V>~ J - ip~ J+2 » 2dz{d z tp)~ J . 

Equation (|5.27j) gives us a system of linear equations to solve at each timestep. The 
coefficients of the system are 3x3 complex matrices because of the gauge fields in the 
equation. Note that although we need the spinor at both timesteps n — 1 and n to find 
ipn+i, the P^ components decouple in such a way that we only need to store the values 
of tp at one timestep, i.e. after using P + tp n _i and P~tp n to find P + tp n+ \ we can then 
forget P + ip n -i and move to the next timestep to solve for P~ip n+ 2- The linear system 
Eq. (|5,27|) is most efficiently solved by LU-decomposition in a way that we will review 
next. 
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Solving the tridiagonal system by LU-decomposition Let us write the linear 
system, Eq. (|5.27|) . as 

MP ± ^ n+1 = £ n , (5.28) 

where £ n as defined by the r.h.s. of Eq. (j5.27j) is known and P ip n +i is unknown. The 
coefficients of the matrix M can be read off Eq. (|5.27|) taking into account the special 
treatment of the boundary in the z-direction. The matrix M is almost tridiagonal, 
meaning that apart from the two exceptional elements mi 3 an d niNN-2 coming from the 
boundary its elements other than — 1), and + 1) are zero. For notational 
simplicity we shall neglect the spinorial (7-matrix) indices; the following calculation 
separates easily for the two components of P^ifj {tp has four complex components, P^ip 
has two that can be separated e.g. into the eigenvectors of 7 .) The LU-decomposition 
of M is defined as follows 



M 



( mu "112 mi 3 

77121 1"22 "123 

m 32 m 33 777,34 



\ 



LU 



( 1 

h 






1 

h 



ITlN-lN-2 mjV-liV-1 tun-in 

m NN _2 


1 



m NN -i m NN ) 
\ 



In In-i 
I d\ U2 u\ 
d 2 u 3 

d 3 U4 



V 



1 / 



d-N-1 UN 

d N J 



(5.29) 



where we have defined N = 2 J + 1. The algorithm for finding the elements of the lower 
and upper triangular matrices L and U is the following 

1. Set d\ = mn, u\ = "ii3> ^2 = "ii2- 

2. Set li = m<2,x/d\ and u 3 = 77123 — hu\. 

3. Set 

Um+l — "lmm+1 for all 777 > 3. 

4. Set d m = m mm - l m -iu m and l m = m m+ i m /d m for all m>2. 



5. The last row is again exceptional: when we have . . . In-2, ■ ■ ■ ^jv-i we set In = 
mNN-2/dN-2, In-i = ("liViV-l — InUn-i) / 'd^-x and d/v = ttiatjv — In-iun- 

Note that this can be done handily in place (without additional memory, replacing the 
elements of M by the elements of L and U). The expressions like 77721 /di actually mean 
"i2i(rfi) _1 ) because d\ is a complex 3 x 3-matrix. 
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Having LU-decomposed the coefficient matrix we can solve the equation by back- 
substitution. The equation we are solving is LUip = £. Let us denote Uip = x\ L% = 
We first find \ by noting that xi = £ii Xn+i = Cn+i — InXm • • • , except for the last 
exceptional element \n = 6v — InXn-2 — In-iXn-i- This is the lower backsubstitu- 
tion. Then follows the upper backsubstitution, starting from the end and progress- 
ing backwards: Vtv = dJ^XN, then ip n = d~ 1 (x n — u n+ iip n+ i) . . . , until the first one 
V'l = d^ixi - u 2 ip 2 ~ uiip 3 ). 

Equation Q5.27JI has now been solved and we can proceed to the next timestep. The 
number of operations in this solution is linearly proportional to the number of lattice 
points, as in an explicit scheme. The constant of proportionality is, however, higher, so 
the algorithm is slower by a constant factor. 
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Chapter 6 

Results 



Let us then summarize the most important findings of Refs. P3I21IS]- 

6.1 Gluon production at central rapidities 

In Ref. m the k^-distribution, the integrated multiplicity and transverse energy of 
gluons produced in the classical field model were calculated. The results corrected an 
error in the normalization in |138[ I193| 1 and were found to agree with the experimental 
data within the limits given by the analysis in Sec. 12.31 Ref. also discussed the 
gauge dependence of the results and the relation to the weak field result of |1311 1139j 
(see Sec. I4.3|) . The phase space density of gluons in this calculation was found to be 
lower than one would expect, meaning that the classical field approximation is really 
applicable only to the very low momentum modes with px < g 2 fJ- Collisions of finite size 
nuclei with different impact parameters and using a modified probability distribution in 
stead of the original McLerran-Venugopalan form Eq. (13.101) were also studied. 



— (1=0.3 GeV, f £ 
li=0.5 GeV, f E 
(1=0.8 GeV, f 




0.3 GeV, f N 
(1=0.5 GeV, f N 
(1=0.8 GeV, f 



Figure 6.1: Left: The dependence of ]e and /at on the field strength parameter g 2 {iR,A- 
Right: The dependence of Je and /at on lattice spacing for different values of [i. 



The energy and multiplicity, parametrized by the two dimensionless ratios 

dE/d V dN/d V 



1 For an explanation of this normalization issue see the erratum |194| . 
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were calculated numerically using the method outlined in Sec. 15.11 These ratios were 
found to be approximately independent of the parameters of the calculation (g, fj,, Ra 
and the lattice spacing a) in the strong field limit (g 2 fi large enough), as can be seen 
from Fig. 16.11 

The values found in Ref. [Q were /e ~ 0.23 and /jv ~ 0.29. These fit in with the 
two phenomenological scenarios outlined in Sec. 12.31 in the following way. Assuming a 
nuclear transverse area of 140 fm 2 , g = 2 and taking g 2 [i = 2.1 GeV one gets 

J )\7"init J pinit 

diV tot „ noo ^^^ 1900 GeVi (62) 



di] drj 

which fits in well with the hydro dynamical expansion scenario outlined in Sec. 12.31 On 
the other hand a choice of g 2 /j. = 1.4 GeV would give 

j Ajinit J cnnit 

tot ~ 500 and —j-- ~ 550 GeV, (6.3) 



drj dr] 

which would be consistent with with the free streaming scenario, assuming a factor of 2 
increase in the multiplicity from the scatterings of the partons and hadronization. 

The CPU-time used for producing the results of Ref. 1. was of the order of a few 
thousand hours on a 2.4 GHz Pentium processor, with comparable amounts of computer 
resources used in the development phase of the code. The calculations were performed 
during the winter 2002-2003 using up to five 600 MHz Alpha EV6 processors in the 
"dynamo" minisupercluster at the Accelerator Laboratory and three desktop 2.4 GHz 
Pentium machines at the Theory Division, both at the University of Helsinki Department 
of Physics. The program consisted of approximately 7000 lines of code in ANSI C in 
addition to general purpose libraries for elementary complex number and SU(3)-matrix 
operations. 



6.2 Rapidity dependence 

The calculations of Ref. were further extended in Ref. |2j ■ The experimentally found 
rapidity (or In 1/x) dependence of the saturation scale discussed in Sec. l3.2l was exploited 
in a simple extension of the calculation of [Q to study the rapidity dependence of gluon 
production. The main assumption of the calculation was the following: in the fully boost 
invariant model described in Sec. 14. II the only (implicit) dependence on rapidity comes 
though the strengths of the classical color sources P(i : 2)- It nas been experimentally 
observed that the saturation scale varies with rapidity as Q 2 (x) ~ x _A . The saturation 
scale in the McLerran-Venugopalan model is proportional, up to a logarithm, to the 
source strength, Q s ~ g 2 //(x lng 2 ^). It was argued in Ref. [2J that for central rapidities 
in heavy ion collisions the dominant x scale that should be used to evaluate the saturation 
scale is given by ~ e ±y for the two nuclei. One can then study gluon production at 
the rapidity y by choosing the source strengths as g 4 /J? ~ e ±Ay . 

The result of Ref. j^j was that the multiplicity of produced gluons can be described by 
very broad Gaussians in y with width a ~ 6. This kind of a Gaussian is much broader 
than the one observed at RHIC energies by the BRAHMS collaboration |204j 2 (see 

2 The rapidity dependence has also been measured by STAR [205) . but the coverage in rapidity is too 
small [\y\ < 0.5) to give a clear interpretation. 
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Figure 6.2: BRAHMS results on charged meason rapidity distributions in 5 % most 
central gold-gold collisions 204 . The dashed lines are Gaussian fits to the data with 
ov+ = 2.25 ± 0.02 (stat) and a„- = 2.29 ± 0.02 (stat). The lower panel shows the mean 
transverse momentum of the mesons for different rapidities. 



Fig. l6,2j) . but coincidentally close to the pQCD+saturation model result |206|I207] 3 . The 
conclusion suggested in Ref. [2] was then that at RHIC energies the rapidity dependence 
of particle production at central rapidities at AuAu-collisions at RHIC energies is not 
dominated by saturation physics. It is perhaps more dependent on the kind of high x 
physics (the (1 — x) 4 -behavior of the gluon structure functions expected from momentum 
sum rules) incorporated e.g. in the "saturation model" calculations of Refs. |171 1931 l9~4"| 
to reproduce the experimentally observed rapidity (or pseudorapidity) dependence. 

The CPU-time used for producing the results of Ref. ;2J was larger than what was 
used for Ref. because the computations had to be done separately for different values 
of the rapidity. The total CPU-time was approximately 5000 hours on a 2.4 GHz Pentium 
processor. The calculations were performed during the late summer of 2004 mostly using 
three desktop 2.4 GHz Pentium machines at the Theory Division of the University of 
Helsinki Department of Physics. 



6.3 Quark pair production: 1+1-dimensional toy model 

In Ref. jS] the calculation of pair production by classical fields was formulated based 
on the theory detailed in Ref. |189j and closely following the calculation in the Abelian 
theory performed in Ref. 188 . It was then argued that the practical numerical solution 

3 Note that according to Ref. 208 hydrodynamical evolution does not change the rapidity distribution 
enough to explain this large a difference between the initial and final states. 
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Figure 6.3: Diagrams contributing to the quark pair prodction amplitude in the 1+1- 
dimensional toy model. 

of the Dirac equation in the forward light cone with the initial condition given at r = 
is most practically performed using as coodinates r and z. The implicit discretization 
scheme required for this numerical solution was then constructed and demonstrated in 
a 1+1-dimensional toy model. 

In the gauge A T = that was used in calculating the classical background field the 
Dirac equation of the 1+1-dimensional toy model is 



d T ijj 



a/t 2 + Z 2 + 7°7 3 2! 



A, 



-7 7 3 <9 z t/> - imefryV) - H°1 3 9—'4' 



(6.4) 



The transverse coordinate dependence has been neglected and the mass of the 1+1 

,2 



k T 2 + 



dimensional theory corresponds to the transverse mass of the full theory, m 2 ff 
m 2 . The 1+1-dimensional toy model in the A T = 0-gauge only contains one component 
of the gauge field, A v . In Ref. |5] this model was studied using two forms for the time 
dependence of this background field, an exponential decay and a Bessel function: 



An = cQ s t.Ji(Q s t). 



(6.5) 



This is the correct time dependence of the perturbative solution to the gauge field 
equations of motion (131 j . 



Amplitude for different values of Q 

Weak fields: cQ = 0.05 m s 

10 



Q s = 1.95 m 

. . . . q[ = 2.0 m 
- - - q[ = 2.05 m 
-- q[ = 2.1 m 
. - Q S = 2.2m 




Amplitude for different values of Q 

Strong fields: cQ s = m s 




Figure 6.4: Absolute value of the quark pair production amplitude for different values 
of the oscillation scale Q s . Left: weak fields, the peaks are at the location given by the 
lowest order perturbative result. Right: strong fields with the same values of Q s - Peaks 
near "threshold" Q s = 2m are shifted. 
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For weak fields c <C 1 the amplitude for quark pair production can be computed 
using the first order in perturbation theory (diagram (a) in Fig. 16, 3|) . The result for the 
oscillating field of Eq. 1)6. 5 jl is a peak at 

2k + k~ = (p + qf = 2m 2 eS (1 + cosh(Ay)) = Q 2 . (6.6) 

As is demonstrated in Fig. 16.41 for weak fields the numerical computation of Ref. |Hj 
reproduces this perturbative peak. For stronger fields the numerical solution sums over 
all the diagrams in Fig. 16.31 and the position of the peak is shifted. The appearence 
of this peak is not a physical effect. The time dependent field of Eq. H6.5|) represents 
off-shell gluons with invariant mass Q 2 decaying into quark pairs, and the amplitude has 
a peak where the invariant mass of the pair equals that of the gluon. In the full 3+1- 
dimensional case the gluons are on shell, and this peak is washed out by integration over 
the relative transverse momentum of the pair. The effect serves here as a verification 
that the numerical method for solving the Dirac equation and projecting out the positive 
energy solutions works. 
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Chapter 7 

Review and outlook 



The classical field approximation is a theoretically appealing way to study QCD as it 
manifests itself in the small x proton or nucleus wavefunction and in relativistic heavy 
ion collisions. The aim of this thesis has been to review how a classical field model can 
be used to understand the initial stages of a heavy ion collision towards thermalization, 
starting from a model of the nuclear wavefunction. It has been seen that the classical 
field model gives a plausible picture of gluon production in the central rapidity region of a 
heavy ion collision. One can numerically calculate the gluon multiplicity and transverse 
energy created in the initial stages of the collision. This calculation was the subject of the 
first article of this thesis pQ, and it was shown that the values obtained are compatible 
with experimental observations within the Bjorken scenario of boost invariant ideal 
hydrodynamical expansion. 

Due to the uncertainty of the appropriate value of the source charge density the 
results of the classical model do not, however, rule out a free-streaming scenario without 
the decrease of energy density due to p dl/-work in the hydronynamical evolution. This 
uncertainty can be reduced in two ways. Firstly observations of the hadron suppression 
at large px, the disappearence of away-side jet-like correlations and elliptic flow favor 
the interpretation in terms of the Bjorken scenario, although hydrodynamical models 
do have difficulties with HBT radii and there have been attempts to explain elliptic flow 
without referring to hydrodynamics |2L)9U137j . Secondly one should in principle be able 
to independently determine the value of the source strength parameter using data from 
deep inelastic scattering. Doing this accurately enough to clearly distinguish between 
g 2 \x = 1.4 GeV and g 2 \x = 2.1 GeV is, however, not easy. 

Inclusion of the longitudinal dimension is important for fully understanding two 
seemingly unrelated issues. One is the thermalization, or specifically the isotropisation, 
of the initial gluonic system. The essential assumption of the hydrodynamical model is 
that the thermal system is isotropic, i.e. has not only transverse but also longitudinal 
pressure. Explaining the creation of longitudinal pressure, if it can be done in the 
classical field model, must necessarily involve a full 3+1-dimensional solution of the 
gauge field equations of motion. The second issue that requires a better treatment of 
the longitudinal dimension is including the effects of JIMWLK evolution. The approach 
used in the second article of the present thesis [2j , consists of varying the initial conditions 
of the same 2+1-dimensional calculation to account for the change in the saturation scale 
at different rapidities. There exists a large body of both analytical and numerical work 
on this equation, which is not yet fully included in this simple approach. 

The question of theoretically understanding of the thermal and chemical equilibra- 
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tion that the experimental results seem to indicate is not as well understood. In the 
third article of this thesis we have formulated a way to calculate the number of quark 
pairs produced in the classical field model. The calculation essentially consists of solving 
the Dirac equation in the background given by the classical gauge fields. Many technical 
issues concerning, in particular, the longitudinal dimension have been solved in a 1+1- 
dimensional toy model. The full 3+1-dimensional case, requiring an extensive numerical 
calculation, has been formulated but not yet implemented numerically. 
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Appendix A 

Notations 



A.l Spacetime 

The metric in the original {t, z, x^j-basis is g^ v = diag(l, — 1, — 1, — 1). Thus did* = 
-didi = -V|. We define the light cone coordinates and proper time and spacetime 
rapidity T,rj as: 

x ± = -^=(t±z) = -^re ±7? 

y/2 y/2 



T = y/t 2 - Z 2 = V2X+X- (A.l) 

1 t + z 1 , x + 

77 = — In = — In — , 

' 2 t-z 2 x- 

which gives the original coordinates in terms of the new ones as 

t = —=(x + + x~) = t cosh 1] 

f (A.2) 

z = —=(x + —x ) = rsinhr?. 
y/2 

The metric in the light cone coordinates is ds 2 = 2dx + x~ — dx^ and the r, rj coordi- 
nates ds 2 = dr 2 — t 2 drj 2 — dx™, giving the invariant integration measure y/[g[d 4 x = 
dzd£d 2 xj* = dx + dx~ d 2 xy = r drj dr d 2 x^. Any vector, in particular the gauge field 
A„, transforms under coordinate transformations in the familiar way 



# = T^" \' = HZ*Ar ( A - 3 ) 



dx^ dx^ 
g x n ^ q x h> 

Straightforward application of Eq. ()A.3|) gives the r, 77-components of a vector as 

A T = A T = (x + A- + X -A + )/t 

_u i (A. 4) 

A.2 Chromodynamics 

We denote the generators of the fundamental representation of the gauge group by t a . 
They are related to the Gell-Mann matrices A a and the Pauli matrices a a by t a = X a /2 
for N c = 3 and t a = cr a /2 for N c = 2. The generators are normalized as 

Trt a t b = -8 ah (A.5) 
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[t a ,t b ] = if abate- (A.6) 



The structure constants are completely antisymmetric in their indices |210| . The color 
indices a, b, c, . . . range from 1 to N c 2 — 1, the dimension of the group. They are written 
here indifferently as superscripts or subscripts. The adjoint representation is generated 
by (T a ) bc = -ifabc- 

TrT a T b = N c 5 ab , or T a T a = iV c l (JVc 2_ 1 ), (A.7) 

where l^ Nc 2_^ is the (N c 2 — 1) x (N c 2 — 1) identity matrix. The adjoint representation 
of an SU(-/V" c )-matrix U = expix a t a is denoted by D^(U) = expix a T a , and is related 
to the fundamental representation by 

2Tr(t a Uh b U) = (DadjCEOU- (A.8) 
We use the following notation for the dimensions and the Casimirs of the representations: 



d F = N c d A = N c z -l 

The action of QCD is 



(A.9) 



S= f d^x^MP-m^ — TrF^F^ (A.10) 
/ 

with the covariant derivative 

D» = & 1 + igA», {All) 

giving the fields trength tensor 

F ^ = -l[D> i ,D v ] =d"A u -d u A» +ig[A' t ,A v ] (A.12) 



or in components 



Fr = d»A v a - WA» - gf abc A^ c . (A.13) 



QCD is invariant under gauge transformations, parametrized by a spacetime depen- 
dent matrix in the gauge group, U(x): 

A» -► UA^rf - -Ud^Ul (A. 14) 

9 

Under these transformations the covariant derivative and the field strength tensor trans- 
form as 

D» -► UD^ and F^ -> UF^rf. (A.15) 
A field configuration that is gauge equivalent to the vacuum is called a pure gauge: 

A* = --Ud^Ul (A.16) 
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Ufa 



Figure A.l: The link matrix U ll (x) connects the lattice site x to x +e M (left). The 
Hermitian conjugate Ul(x) connects lattice points in the negative direction, going from 
x +e„ to x. 

A. 3 Gauge fields on the lattice 

The fundamental object of lattice gauge field theory is the link matrix Uaix) connecting 
the point x to the neighboring lattice site x+e^ (see Fig. lA.ljl . The Hermitian conjugate 
of the link matrix reverses the direction, connecting x+e^ to x. In terms of the continuum 
gauge field we identify the link matrix as Uu(x) = e l9aA ^^ 

A gauge transformation on the lattice transforms the link matrices as 

U^x) -> y(x)C/ M (x)W(x+e M ), (A.17) 

and the fermion fields as 

ip(x) -> V(x)i/>(x), 4>(x) -> $(x)V\ (A.18) 

leaving the following combination gauge-invariant 

^{x^U^x^x+e^). (A.19) 

A pure gauge field is one that can be gauge transformed to zero, meaning that the link 
matrices can be transformed to the indentity matrix: 

U^x) = V(x)V\x+eJ. (A.20) 

A gauge invariant action can be defined using the plaquette 

U^{x) = U^x)U v (x+ejUl(x+e v )Ut(x). (A.21) 

With a subscript — fi we denote the negative //-direction. Thus a plaquette around the 
same contour, but based on a different point (see Fig. IA.2|) can be written e.g. as 

U v ^{x = U v {x+e^Ul(x+e v )Ul(x)U„(x). (A.22) 

As can be seen from its definition and the transformation properties of the link matrices, 
the plaquette gauge transforms at its base point 

U^{x) - V(x)U^(x)V\x). (A.23) 

Because of cyclicity, the trace of the plaquette is gauge invariant and independent of the 
base point 

TV U^ v {x) = Tr U v ^(x +e^) = Tr U^ v {x +e^+e u ) = Tr U.^(x +e„). (A.24) 
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x+e v x+e fl +e v x+e v x+e^ + e^ 




U v (x +e M ) 




U v (x +e M ) 



x U^{x) x+e^ x Unix) x+e^ 



Figure A. 2: The plaquettes U^.^(x) (left) and [/^-^(x +e M ) (right) circle the same 
contour, but gauge transform on different points, x and x +e M respectively. 

In the continuum limit the plaquette becomes 

and the pure gauge action S = — | Tr F^ u can be represented by the Wilson action |211j 

^-^EEf 1 -^ 11 ^)- (A - 26) 

In Minkowski space one must choose a different sign for the 0, z-terms to reproduce the 
right metric. 
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